#region Copyright Notice
// 
// 
//  
//  Copyright (c) 2009-2011 Anindya Chatterjee
//  
//  Licensed under the Apache License, Version 2.0 (the "License");
//  you may not use this file except in compliance with the License.
//  You may obtain a copy of the License at
//  
//      http://www.apache.org/licenses/LICENSE-2.0
//  
//  Unless required by applicable law or agreed to in writing, software
//  distributed under the License is distributed on an "AS IS" BASIS,
//  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
//  See the License for the specific language governing permissions and
//  limitations under the License.
// 
// 
#endregion
using System;
using System.Collections.Generic;
using AbstractClass.Extensions;
using AbstractClass.Linq;

namespace AbstractClass
{
    /// <summary>
    /// Provides a set of mathematical functions which are not present in <see cref="Math"/> class.
    /// </summary>
    public static class QMath
    {
        private const int MinLength = 2;
        private const int MaxLength = 16384;
        private const int MinBits = 1;
        private const int MaxBits = 14;
        private static readonly int[][] ReversedBits = new int[MaxBits][];
        private static readonly ComplexNumber[,][] ComplexRotation = new ComplexNumber[MaxBits,2][];

        /// <summary>
        /// Calculates the inverse sine hyperbolic of a specified number.
        /// </summary>
        /// <param name="number">The number.</param>
        /// <returns>Inverse sine hyperbolic of the specified number.</returns>
        public static double Asinh(double number)
        {
            return Math.Log(number + Math.Sqrt(number.ToThePower(2) + 1));
        }

        /// <summary>
        /// Calculates the inverse cos hyperbolic of a specified number.
        /// </summary>
        /// <param name="number">The number.</param>
        /// <returns>Inverse cos hyperbolic of the specified number.</returns>
        public static double Acosh(double number)
        {
            return Math.Log(number + Math.Sqrt(number.ToThePower(2) - 1));
        }

        /// <summary>
        /// Calculates the inverse tan hyperbolic of a specified number.
        /// </summary>
        /// <param name="number">The number.</param>
        /// <returns>Inverse tan hyperbolic of the specified number.</returns>
        public static double Atanh(double number)
        {
            return (Math.Log(1 + number) - Math.Log(1 - number)) / 2;
        }


        /// <summary>
        /// Calculates the Error function of a specified number.
        /// </summary>
        /// <param name="number">The number.</param>
        /// <returns>The value of the error function of the specified number.</returns>
        public static double Erf(double number)
        {
            const double a1 = 0.254829592;
            const double a2 = -0.284496736;
            const double a3 = 1.421413741;
            const double a4 = -1.453152027;
            const double a5 = 1.061405429;
            const double p = 0.3275911;

            int sign = 1;
            if (number < 0)
                sign = -1;
            number = Math.Abs(number);

            double t = 1.0 / (1.0 + p * number);
            double y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * Math.Exp(-number * number);

            return sign * y;
        }

        /// <summary>
        /// Calculates the cumulative density function (CDF) of a standard normal (Gaussian) random variable.
        /// </summary>
        /// <param name="number">The standard normal (Gaussian) random variable.</param>
        /// <returns>Cumulative density function (CDF) of the standard normal (Gaussian) random variable.</returns>
        public static double Phi(double number)
        {
            const double a1 = 0.254829592;
            const double a2 = -0.284496736;
            const double a3 = 1.421413741;
            const double a4 = -1.453152027;
            const double a5 = 1.061405429;
            const double p = 0.3275911;

            int sign = 1;
            if (number < 0)
                sign = -1;
            number = Math.Abs(number) / Math.Sqrt(2.0);

            double t = 1.0 / (1.0 + p * number);
            double y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * Math.Exp(-number * number);

            return 0.5 * (1.0 + sign * y);
        }

        /// <summary>
        /// Calculates the cumulative density function (CDF) inverse of a standard normal (Gaussian) random variable.
        /// </summary>
        /// <param name="number">The standard normal (Gaussian) random variable.</param>
        /// <returns>Cumulative density function (CDF) inverse of the standard normal (Gaussian) random variable.</returns>
        /// <exception cref="ArgumentOutOfRangeException">Specified number is less than equal to 0 or greater than or equal to 1.</exception>
        public static double PhiInverse(double number)
        {
            if (number <= 0.0 || number >= 1.0)
            {
                string msg = String.Format(LanguageResource.QMathPhiInverseInvalidInput, number);
                throw new ArgumentOutOfRangeException(msg);
            }

            return number < 0.5
                       ? -RationalApproximation(Math.Sqrt(-2.0 * Math.Log(number)))
                       : RationalApproximation(Math.Sqrt(-2.0 * Math.Log(1.0 - number)));
        }

        /// <summary>
        /// Calculates the rational approximation of the <paramref name="number"/>.
        /// </summary>
        /// <param name="number">The number.</param>
        /// <returns>Rational approximation of the <paramref name="number"/>.</returns>
        public static double RationalApproximation(double number)
        {
            // Abramowitz and Stegun formula 26.2.23.
            // The absolute value of the error should be less than 4.5 e-4.
            double[] c = {2.515517, 0.802853, 0.010328};
            double[] d = {1.432788, 0.189269, 0.001308};
            return number - ((c[2] * number + c[1]) * number + c[0]) /
                   (((d[2] * number + d[1]) * number + d[0]) * number + 1.0);
        }

        /// <summary>
        /// Computes log(1 + <paramref name="number"/>) without losing precision for small values of <paramref name="number"/>.
        /// </summary>
        /// <param name="number">The number.</param>
        /// <returns>The value of log(1 + <paramref name="number"/>).</returns>
        /// <exception cref="ArgumentOutOfRangeException">The specified number is less than or equal to -1.</exception>
        public static double LogOnePlusX(double number)
        {
            if (number <= -1.0)
            {
                string msg = String.Format(LanguageResource.QMathLogOnePlusXInvalidInput, number);
                throw new ArgumentOutOfRangeException(msg);
            }

            if (Math.Abs(number) > 1e-4)
            {
                return Math.Log(1.0 + number);
            }

            // Use Taylor approx. log(1 + x) = x - x^2/2 with error roughly x^3/3
            // Since |x| < 10^-4, |x|^3 < 10^-12, relative error less than 10^-8
            return (-0.5 * number + 1.0) * number;
        }


        /// <summary>
        /// Calculates exp(<paramref name="number"/>) - 1 without losing precision for small values of <paramref name="number"/>.
        /// </summary>
        /// <param name="number">The number.</param>
        /// <returns>The value of exp(<paramref name="number"/>) - 1.</returns>
        public static double ExpMinusOne(double number)
        {
            return Math.Abs(number) < 1e-5 ? number + 0.5 * number * number : Math.Exp(number) - 1.0;
        }

        /// <summary>
        /// Determines whether a specified number is prime.
        /// </summary>
        /// <param name="number">The number.</param>
        /// <returns>
        /// 	<c>true</c> if the specified number is prime; otherwise, <c>false</c>.
        /// </returns>
        public static bool IsPrime(int number)
        {
            switch ((number & 1))
            {
                case 0:
                    return number == 2;
                default:
                    double sqrt = Math.Sqrt(number);
                    for (int i = 3; i <= sqrt; i += 2)
                    {
                        if ((number % i) == 0)
                        {
                            return false;
                        }
                    }
                    return true;
            }
        }

        /// <summary>
        /// Calculates Discrete Fourier Transform of the specified data.
        /// </summary>
        /// <param name="data">The collection of <see cref="ComplexNumber"/> data.</param>
        /// <param name="direction">The Fourier transformation direction.</param>
        /// <returns>An array of <see cref="ComplexNumber"/> containing the transformed values.</returns>
        public static ComplexNumber[] DiscreteFourierTransform(ComplexNumber[] data, FourierDirection direction)
        {
            int n = data.Length;
            var dst = new ComplexNumber[n];

            for (int i = 0; i < n; i++)
            {
                dst[i] = ComplexNumber.Zero;
                double arg = -(int) direction * 2.0 * Math.PI * i / n;

                for (int j = 0; j < n; j++)
                {
                    double cos = Math.Cos(j * arg);
                    double sin = Math.Sin(j * arg);

                    dst[i].Real += (float) (data[j].Real * cos - data[j].Imaginary * sin);
                    dst[i].Imaginary += (float) (data[j].Real * sin + data[j].Imaginary * cos);
                }
            }

            if (direction == FourierDirection.Forward)
            {
                for (int i = 0; i < n; i++)
                {
                    dst[i].Real = dst[i].Real / n;
                    dst[i].Imaginary = dst[i].Imaginary / n;
                }
            }
            else
            {
                for (int i = 0; i < n; i++)
                {
                    dst[i].Real = dst[i].Real;
                    dst[i].Imaginary = dst[i].Imaginary;
                }
            }
            return dst;
        }


        /// <summary>
        /// Calculates 2-dimensional Discrete Fourier Transform of the specified data.
        /// </summary>
        /// <param name="data">The collection of <see cref="ComplexNumber"/> data.</param>
        /// <param name="direction">The Fourier transformation direction.</param>
        /// <returns>An array of <see cref="ComplexNumber"/> containing the transformed values.</returns>
        [CLSCompliant(false)]
        public static ComplexNumber[,] DiscreteFourierTransform(ComplexNumber[,] data, FourierDirection direction)
        {
            int n = data.GetLength(0);
            int m = data.GetLength(1);
            double arg, cos, sin;
            var dst = new ComplexNumber[Math.Max(n, m)];
            var dummyData = new ComplexNumber[n,m];
            data.CopyTo(dummyData);

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < m; j++)
                {
                    dst[j] = ComplexNumber.Zero;

                    arg = -(int) direction * 2.0 * Math.PI * j / m;

                    for (int k = 0; k < m; k++)
                    {
                        cos = Math.Cos(k * arg);
                        sin = Math.Sin(k * arg);

                        dst[j].Real += (float) (dummyData[i, k].Real * cos - dummyData[i, k].Imaginary * sin);
                        dst[j].Imaginary += (float) (dummyData[i, k].Real * sin + dummyData[i, k].Imaginary * cos);
                    }
                }

                if (direction == FourierDirection.Forward)
                {
                    for (int j = 0; j < m; j++)
                    {
                        dummyData[i, j].Real = dst[j].Real / m;
                        dummyData[i, j].Imaginary = dst[j].Imaginary / m;
                    }
                }
                else
                {
                    for (int j = 0; j < m; j++)
                    {
                        dummyData[i, j].Real = dst[j].Real;
                        dummyData[i, j].Imaginary = dst[j].Imaginary;
                    }
                }
            }

            for (int j = 0; j < m; j++)
            {
                for (int i = 0; i < n; i++)
                {
                    dst[i] = ComplexNumber.Zero;

                    arg = -(int) direction * 2.0 * Math.PI * i / n;

                    for (int k = 0; k < n; k++)
                    {
                        cos = Math.Cos(k * arg);
                        sin = Math.Sin(k * arg);

                        dst[i].Real += (float) (dummyData[k, j].Real * cos - dummyData[k, j].Imaginary * sin);
                        dst[i].Imaginary += (float) (dummyData[k, j].Real * sin + dummyData[k, j].Imaginary * cos);
                    }
                }

                if (direction == FourierDirection.Forward)
                {
                    for (int i = 0; i < n; i++)
                    {
                        dummyData[i, j].Real = dst[i].Real / n;
                        dummyData[i, j].Imaginary = dst[i].Imaginary / n;
                    }
                }
                else
                {
                    for (int i = 0; i < n; i++)
                    {
                        dummyData[i, j].Real = dst[i].Real;
                        dummyData[i, j].Imaginary = dst[i].Imaginary;
                    }
                }
            }
            return dummyData;
        }


        /// <summary>
        /// Calculates Fast Fourier Transform of the specified data.
        /// </summary>
        /// <param name="data">The collection of <see cref="ComplexNumber"/> data.</param>
        /// <param name="direction">The Fourier transformation direction.</param>
        /// <returns>An array of <see cref="ComplexNumber"/> containing the transformed values.</returns>
        public static ComplexNumber[] FastFourierTransform(ComplexNumber[] data, FourierDirection direction)
        {
            int n = data.Length;
            var m = (int) (Math.Log(n, 2));
            var dummyData = new ComplexNumber[n];
            data.CopyTo(dummyData, 0);

            ReorderData(dummyData);

            int tn = 1;

            for (int k = 1; k <= m; k++)
            {
                ComplexNumber[] rotation = GetComplexRotation(k, direction);

                int tm = tn;
                tn <<= 1;

                for (int i = 0; i < tm; i++)
                {
                    ComplexNumber t = rotation[i];

                    for (int even = i; even < n; even += tn)
                    {
                        int odd = even + tm;
                        ComplexNumber ce = dummyData[even];
                        ComplexNumber co = dummyData[odd];

                        var tr = (float) (co.Real * t.Real - co.Imaginary * t.Imaginary);
                        var ti = (float) (co.Real * t.Imaginary + co.Imaginary * t.Real);

                        dummyData[even].Real += tr;
                        dummyData[even].Imaginary += ti;

                        dummyData[odd].Real = ce.Real - tr;
                        dummyData[odd].Imaginary = ce.Imaginary - ti;
                    }
                }
            }

            if (direction == FourierDirection.Forward)
            {
                for (int i = 0; i < n; i++)
                {
                    dummyData[i].Real /= n;
                    dummyData[i].Imaginary /= n;
                }
            }

            return dummyData;
        }


        /// <summary>
        /// Calculates 2-dimensional Fast Fourier Transform of the specified data.
        /// </summary>
        /// <param name="data">The collection of <see cref="ComplexNumber"/> data.</param>
        /// <param name="direction">The Fourier transformation direction.</param>
        /// <returns>An array of <see cref="ComplexNumber"/> containing the transformed values.</returns>
        /// <exception cref="ArgumentException">The number of elements in any dimension of the array satisfies any of the following conditions
        /// <list type="">
        /// <item>Less than 2</item>
        /// <item>Can not be represented as a power of 2</item>
        /// <item>Greater than 16384</item>
        /// </list>
        /// </exception>
        [CLSCompliant(false)]
        public static ComplexNumber[,] FastFourierTransform(ComplexNumber[,] data, FourierDirection direction)
        {
            int k = data.GetLength(0);
            int n = data.GetLength(1);
            var dummyData = new ComplexNumber[k,n];
            data.CopyTo(dummyData);

            if ((!IsPowerOf2(k)) || (!IsPowerOf2(n)) ||
                (k < MinLength) || (k > MaxLength) ||
                (n < MinLength) || (n > MaxLength))
            {
                throw new ArgumentException();
            }

            var row = new ComplexNumber[n];

            for (int i = 0; i < k; i++)
            {
                for (int j = 0; j < n; j++)
                    row[j] = dummyData[i, j];
                FastFourierTransform(row, direction);
                for (int j = 0; j < n; j++)
                    dummyData[i, j] = row[j];
            }

            var col = new ComplexNumber[k];

            for (int j = 0; j < n; j++)
            {
                for (int i = 0; i < k; i++)
                    col[i] = dummyData[i, j];
                FastFourierTransform(col, direction);
                for (int i = 0; i < k; i++)
                    dummyData[i, j] = col[i];
            }

            return dummyData;
        }

        private static int[] GetReversedBits(int numberOfBits)
        {
            if ((numberOfBits < MinBits) || (numberOfBits > MaxBits))
                throw new ArgumentOutOfRangeException();

            if (ReversedBits[numberOfBits - 1] == null)
            {
                int n = Pow2(numberOfBits);
                var rBits = new int[n];

                for (int i = 0; i < n; i++)
                {
                    int oldBits = i;
                    int newBits = 0;

                    for (int j = 0; j < numberOfBits; j++)
                    {
                        newBits = (newBits << 1) | (oldBits & 1);
                        oldBits = (oldBits >> 1);
                    }
                    rBits[i] = newBits;
                }
                ReversedBits[numberOfBits - 1] = rBits;
            }
            return ReversedBits[numberOfBits - 1];
        }


        private static ComplexNumber[] GetComplexRotation(int numberOfBits, FourierDirection direction)
        {
            int directionIndex = (direction == FourierDirection.Forward) ? 0 : 1;

            if (ComplexRotation[numberOfBits - 1, directionIndex] == null)
            {
                int n = 1 << (numberOfBits - 1);
                float uR = 1.0f;
                float uI = 0.0f;
                double angle = Math.PI / n * (int) direction;
                var wR = (float) Math.Cos(angle);
                var wI = (float) Math.Sin(angle);
                var rotation = new ComplexNumber[n];

                for (int i = 0; i < n; i++)
                {
                    rotation[i] = new ComplexNumber(uR, uI);
                    float t = uR * wI + uI * wR;
                    uR = uR * wR - uI * wI;
                    uI = t;
                }

                ComplexRotation[numberOfBits - 1, directionIndex] = rotation;
            }
            return ComplexRotation[numberOfBits - 1, directionIndex];
        }

        private static void ReorderData(IList<ComplexNumber> data)
        {
            int len = data.Count;

            if ((len < MinLength) || (len > MaxLength) || (!IsPowerOf2(len)))
                throw new ArgumentException();

            int[] rBits = GetReversedBits((int) (Math.Ceiling(Math.Log(len, 2))));

            for (int i = 0; i < len; i++)
            {
                int s = rBits[i];

                if (s > i)
                {
                    ComplexNumber t = data[i];
                    data[i] = data[s];
                    data[s] = t;
                }
            }
        }

        private static bool IsPowerOf2(int x)
        {
            return (x & (x - 1)) == 0;
        }

        private static int Pow2(int exp)
        {
            return ((exp >= 0) && (exp <= 30)) ? (1 << exp) : 0;
        }

        /// <summary>
        /// Calculates the square root of a specified real number.
        /// </summary>
        /// <param name="number">The real number.</param>
        /// <returns>
        /// A <see cref="ComplexNumber"/> representing the square root of a real number.
        /// </returns>
        public static ComplexNumber Sqrt(double number)
        {
            return number < 0 ? new ComplexNumber(0, Math.Sqrt(number * -1)) : new ComplexNumber(Math.Sqrt(number));
        }

        /// <summary>
        /// Calculates the Cos of the specified <see cref="ComplexNumber"/>.
        /// </summary>
        /// <param name="z">The <see cref="ComplexNumber"/>.</param>
        /// <returns>Cos of the specified <see cref="ComplexNumber"/>.</returns>
        public static ComplexNumber Cos(ComplexNumber z)
        {
            ComplexNumber z1 = Exp(new ComplexNumber(-z.Imaginary, z.Real));
            ComplexNumber z2 = Exp(new ComplexNumber(z.Imaginary, -z.Real));

            return new ComplexNumber((0.5 * (z1.Real + z2.Real)), (0.5 * (z1.Imaginary + z2.Imaginary)));
        }

        /// <summary>
        /// Calculates the Cosh of the specified <see cref="ComplexNumber"/>.
        /// </summary>
        /// <param name="z">The <see cref="ComplexNumber"/>.</param>
        /// <returns>Cosh of the specified <see cref="ComplexNumber"/>.</returns>
        public static ComplexNumber Cosh(ComplexNumber z)
        {
            ComplexNumber z1 = Exp(z);
            ComplexNumber z2 = Exp(new ComplexNumber(-z.Real, -z.Imaginary));

            return new ComplexNumber((0.5 * (z1.Real + z2.Real)), (0.5 * (z1.Imaginary + z2.Imaginary)));
        }

        /// <summary>
        /// Calculates the Sin of the specified <see cref="ComplexNumber"/>.
        /// </summary>
        /// <param name="z">The <see cref="ComplexNumber"/>.</param>
        /// <returns>Sin of the specified <see cref="ComplexNumber"/>.</returns>
        public static ComplexNumber Sin(ComplexNumber z)
        {
            ComplexNumber z1 = Exp(new ComplexNumber(-z.Imaginary, z.Real));
            ComplexNumber z2 = Exp(new ComplexNumber(z.Imaginary, -z.Real));

            return new ComplexNumber((0.5 * (z1.Imaginary - z2.Imaginary)), (0.5 * (z2.Real - z1.Real)));
        }

        /// <summary>
        /// Calculates the Sinh of the specified <see cref="ComplexNumber"/>.
        /// </summary>
        /// <param name="z">The <see cref="ComplexNumber"/>.</param>
        /// <returns>Sinh of the specified <see cref="ComplexNumber"/>.</returns>
        public static ComplexNumber Sinh(ComplexNumber z)
        {
            ComplexNumber z1 = Exp(z);
            ComplexNumber z2 = Exp(new ComplexNumber(-z.Real, -z.Imaginary));

            return new ComplexNumber((0.5 * (z1.Real - z2.Real)), (0.5 * (z1.Imaginary - z2.Imaginary)));
        }

        /// <summary>
        /// Calculates the Tan of the specified <see cref="ComplexNumber"/>.
        /// </summary>
        /// <param name="z">The <see cref="ComplexNumber"/>.</param>
        /// <returns>Tan of the specified <see cref="ComplexNumber"/>.</returns>
        public static ComplexNumber Tan(ComplexNumber z)
        {
            return Sin(z) / Cos(z);
        }

        /// <summary>
        /// Calculates the Tanh of the specified <see cref="ComplexNumber"/>.
        /// </summary>
        /// <param name="z">The <see cref="ComplexNumber"/>.</param>
        /// <returns>Tanh of the specified <see cref="ComplexNumber"/>.</returns>
        public static ComplexNumber Tanh(ComplexNumber z)
        {
            return Sinh(z) / Cosh(z);
        }

        /// <summary>
        /// Calculates the Cot of the specified <see cref="ComplexNumber"/>.
        /// </summary>
        /// <param name="z">The <see cref="ComplexNumber"/>.</param>
        /// <returns>Cot of the specified <see cref="ComplexNumber"/>.</returns>
        public static ComplexNumber Cot(ComplexNumber z)
        {
            return Cos(z) / Sin(z);
        }

        /// <summary>
        /// Calculates the Coth of the specified <see cref="ComplexNumber"/>.
        /// </summary>
        /// <param name="z">The <see cref="ComplexNumber"/>.</param>
        /// <returns>Coth of the specified <see cref="ComplexNumber"/>.</returns>
        public static ComplexNumber Coth(ComplexNumber z)
        {
            return (Exp(2 * z) + 1) / (Exp(2 * z) - 1);
        }

        /// <summary>
        /// Calculates the Sech of the specified <see cref="ComplexNumber"/>.
        /// </summary>
        /// <param name="z">The <see cref="ComplexNumber"/>.</param>
        /// <returns>Sech of the specified <see cref="ComplexNumber"/>.</returns>
        public static ComplexNumber Sech(ComplexNumber z)
        {
            return Inverse(Cosh(z));
        }

        /// <summary>
        /// Calculates the Cosech of the specified <see cref="ComplexNumber"/>.
        /// </summary>
        /// <param name="z">The <see cref="ComplexNumber"/>.</param>
        /// <returns>Cosech of the specified <see cref="ComplexNumber"/>.</returns>
        public static ComplexNumber Cosech(ComplexNumber z)
        {
            return Inverse(Sinh(z));
        }

        /// <summary>
        /// Calculates the e to the power of the specified <see cref="ComplexNumber"/>.
        /// </summary>
        /// <param name="z">The <see cref="ComplexNumber"/>.</param>
        /// <returns>e to the power of the specified <see cref="ComplexNumber"/>.</returns>
        public static ComplexNumber Exp(ComplexNumber z)
        {
            double value = Math.Exp(z.Real);

            return new ComplexNumber(
                (value * Math.Cos(z.Imaginary)),
                (value * Math.Sin(z.Imaginary)));
        }

        /// <summary>
        /// Calculates the natural logarithm of the specified <see cref="ComplexNumber"/>.
        /// </summary>
        /// <param name="z">The <see cref="ComplexNumber"/>.</param>
        /// <returns>Natural logarithm of the specified <see cref="ComplexNumber"/>.</returns>
        public static ComplexNumber Log(ComplexNumber z)
        {
            return new ComplexNumber(Math.Log(z.Modulus), z.Argument);
        }

        /// <summary>
        /// Calculates the base-10 logarithm of the specified <see cref="ComplexNumber"/>.
        /// </summary>
        /// <param name="z">The <see cref="ComplexNumber"/>.</param>
        /// <returns>Base-10 logarithm of the specified <see cref="ComplexNumber"/>.</returns>
        public static ComplexNumber Log10(ComplexNumber z)
        {
            const double log10 = 2.3025850929940459;

            ComplexNumber value = Log(z);

            value.Real = (value.Real / log10);
            value.Imaginary = (value.Imaginary / log10);

            return value;
        }

        /// <summary>
        /// Calculates the square root of the specified <see cref="ComplexNumber"/>.
        /// </summary>
        /// <param name="complexNumber">The complex number.</param>
        /// <returns>
        /// A <see cref="ComplexNumber"/> representing the square root of the complex number.
        /// </returns>
        public static ComplexNumber Sqrt(ComplexNumber complexNumber)
        {
            return ComplexNumber.Polar(Math.Sqrt(complexNumber.Modulus), (complexNumber.Argument * 0.5));
        }

        /// <summary>
        /// Inverses the specified <see cref="ComplexNumber"/>.
        /// </summary>
        /// <param name="complexNumber">The complex number.</param>
        /// <returns>Inverse of the specified complex number.</returns>
        public static ComplexNumber Inverse(ComplexNumber complexNumber)
        {
            return new ComplexNumber(complexNumber.Real / complexNumber.SquaredMagnitude,
                                     -complexNumber.Imaginary / complexNumber.SquaredMagnitude);
        }

        /// <summary>
        /// Returns the larger of two complex numbers.
        /// </summary>
        /// <param name="a">The first of two complex numbers to compare.</param>
        /// <param name="b">The second of two complex numbers to compare.</param>
        /// <returns>The larger of two complex numbers.</returns>
        public static ComplexNumber Max(ComplexNumber a, ComplexNumber b)
        {
            return a.Modulus > b.Modulus ? a : b;
        }
    }
}